Graphical model experiment

The 8th KIAS CAC Summer School

Machine Learning Practice III - Graphical models example

Author: Yung-Kyun Noh, Ph.D.

We show the application of graphical model for Gaussian data. We consider each class has the casuality graph,

Data generation

\begin{eqnarray} \left(\begin{array}{c} X_1 \\ X_2 \end{array}\right) \sim \mathcal{N}(\mu_{12}, \Sigma_{12}), \quad \left(\begin{array}{c} X_7 \\ X_8 \\ X_9 \\ X_{10} \end{array}\right) \sim \mathcal{N}(\mu_{7890}, \Sigma_{7890}), \end{eqnarray}\begin{eqnarray} \left(\begin{array}{c} X_3 \\ X_4 \\ X_5 \end{array}\right) &=& W_{345,12} \left(\begin{array}{c} X_1 \\ X_2 \end{array}\right) + \gamma\left(\begin{array}{c} \epsilon_3 \\ \epsilon_4 \\ \epsilon_5 \end{array}\right) \\ X_6 &=& W_{6,345} \left(\begin{array}{c} X_3 \\ X_4 \\ X_5 \end{array}\right) + W_{6,7890} \left(\begin{array}{c} X_7 \\ X_8 \\ X_9 \\ X_{10} \end{array}\right) + \gamma \epsilon_6 \end{eqnarray}\begin{eqnarray} \epsilon_3, \epsilon_4, \epsilon_5, \epsilon_6 \sim \mathcal{N}(0, 1) \end{eqnarray}

Parameter estimation

\begin{eqnarray} P(X_1,\ldots,X_{10}) &=& P(X_1,X_2)P(X_3,X_4,X_5|X_1,X_2)P(X_7,X_8,X_9,X_{10})P(X_6|X_3,X_4,X_5,X_7,\ldots,X_{10})\\ &=& P(X_1,X_2)\frac{P(X_1,X_2,X_3,X_4,X_5)}{P(X_1,X_2)}P(X_7,X_8,X_9,X_{10})\frac{P(X_3,\ldots,X_{10})}{P(X_3,X_4,X_5,X_7,\ldots,X_{10})} \\ &=& P(X_1,X_2)\frac{P(X_1,X_2,X_3,X_4,X_5)}{P(X_1,X_2)}P(X_7,X_8,X_9,X_{10})\frac{P(X_3,\ldots,X_{10})}{P(X_3,X_4,X_5)P(X_7,\ldots,X_{10})} \\ &=& P(X_1,X_2,X_3,X_4,X_5)\frac{P(X_3,X_4,X_5,X_6,X_7,X_8,X_9,X_{10})}{P(X_3,X_4,X_5)} \end{eqnarray}
In [1]:
import numpy as np

%pylab inline

Populating the interactive namespace from numpy and matplotlib

Make true density parameters

In [2]:
dim = 10;

# True density
trueMean_1 = np.zeros(dim)
trueMean_2 = np.concatenate(([1],np.zeros(dim-1)))
print trueMean_1
print trueMean_2

# covariance matrix of X1, X2 (root nodes)
curDim = 2
trueCov12_1 = np.random.randn(curDim,curDim)
trueCov12_1 = np.dot(trueCov12_1.T, trueCov12_1)
trueCov12_2 = np.random.randn(curDim,curDim)
trueCov12_2 = trueCov12_1 + .001*np.dot(trueCov12_2.T, trueCov12_2)

print trueCov12_1
print trueCov12_2

# covariance matrix of X7, X8, X9, X10 (root nodes)
curDim = 4;
trueCov7890_1 = np.random.randn(curDim,curDim)
trueCov7890_1 = np.dot(trueCov7890_1.T, trueCov7890_1)
trueCov7890_2 = np.random.randn(curDim,curDim)
trueCov7890_2 = trueCov7890_1 + .1*np.dot(trueCov7890_2.T, trueCov7890_2)

print trueCov7890_1
print trueCov7890_2

[ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
[ 1.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
[[ 0.66990231  1.05063391]
 [ 1.05063391  2.23574801]]
[[ 0.67841322  1.04711639]
 [ 1.04711639  2.2385393 ]]
[[ 2.95261466  0.89918369 -0.71586731  0.00907099]
 [ 0.89918369  2.37577021  0.96099359 -0.29608912]
 [-0.71586731  0.96099359  1.18659954  1.00897774]
 [ 0.00907099 -0.29608912  1.00897774  4.05881863]]
[[ 3.16863259  0.94613626 -0.65645078  0.05682842]
 [ 0.94613626  2.5625525   1.17490649 -0.49851749]
 [-0.65645078  1.17490649  1.48813987  0.73875608]
 [ 0.05682842 -0.49851749  0.73875608  4.36524807]]

Generate data according to the true parameters

In [3]:
# data generation
datanum = 11
tstdatanum = 1000
In [4]:
# root node X1, X2
curDim = 2

X12_1 = np.random.multivariate_normal(trueMean_1[0:2], trueCov12_1, datanum)
tstX12_1 = np.random.multivariate_normal(trueMean_1[0:2], trueCov12_1, tstdatanum)

X12_2 = np.random.multivariate_normal(trueMean_2[0:2], trueCov12_2, datanum)
tstX12_2 = np.random.multivariate_normal(trueMean_2[0:2], trueCov12_2, tstdatanum)

# root node X7, X8, X9, X10
curDim = 4
X7890_1 = np.random.multivariate_normal(trueMean_1[6:10], trueCov7890_1, datanum)
tstX7890_1 = np.random.multivariate_normal(trueMean_1[6:10], trueCov7890_1, tstdatanum)
X7890_2 = np.random.multivariate_normal(trueMean_2[6:10], trueCov7890_2, datanum)
tstX7890_2 = np.random.multivariate_normal(trueMean_2[6:10], trueCov7890_2, tstdatanum)
In [5]:
# Non-root node X3, X4, X5
curDim = 3
curW = np.random.randn(2,curDim)
#print curW
X345_1 = np.dot(X12_1,curW) + .1*np.random.randn(datanum,curDim) # non-isotropic noise can be added
#print X345_1
tstX345_1 = np.dot(tstX12_1,curW) + .1*np.random.randn(tstdatanum,curDim) # non-isotropic noise can be added
curW = curW + .1*np.random.randn(2,curDim)
X345_2 = np.dot(X12_2,curW) + .1*np.random.randn(datanum,curDim)   # non-isotropic noise can be added
tstX345_2 = np.dot(tstX12_2,curW) + .1*np.random.randn(tstdatanum,curDim)  # non-isotropic noise can be added

# Non-root node X6
curDim = 1
curW = np.random.randn(7,curDim)
#print np.hstack([X345_1,X7890_1])
X6_1 = np.dot(np.hstack([X345_1,X7890_1]),curW) + .1*np.random.randn(datanum,curDim) # non-isotropic noise can be added
tstX6_1 = np.dot(np.hstack([tstX345_1,tstX7890_1]),curW) + .1*np.random.randn(tstdatanum,curDim) # non-isotropic noise can be added
curW = curW + .1*np.random.randn(7,curDim)
X6_2 = np.dot(np.hstack([X345_2,X7890_2]),curW) + .1*np.random.randn(datanum,curDim) # non-isotropic noise can be added
tstX6_2 = np.dot(np.hstack([tstX345_2,tstX7890_2]),curW) + .1*np.random.randn(tstdatanum,curDim) # non-isotropic noise can be added
In [6]:
X1 = np.hstack([X12_1,X345_1,X6_1,X7890_1])
X2 = np.hstack([X12_2,X345_2,X6_2,X7890_2])
tstX1 = np.hstack([tstX12_1,tstX345_1,tstX6_1,tstX7890_1])
tstX2 = np.hstack([tstX12_2,tstX345_2,tstX6_2,tstX7890_2])
In [7]:
# Classification with joint density function
# X1, X2 \in R^{dim x datanum}
estM1 = sum(X1,axis=0)/datanum
estM2 = sum(X2,axis=0)/datanum

estCov1 = np.dot(X1.T,X1)/datanum - np.outer(estM1, estM1)
estCov2 = np.dot(X2.T,X2)/datanum - np.outer(estM2, estM2)
In [8]:
def getLogGaussian(X, mu, Sig):
    ''' evalute log probability of Gaussian distribution
    given the Gaussian distribution model (mean and covariance).'''
    # X: datanum x dim
    nData, nDim = X.shape
    detSig = np.linalg.det(Sig)
    invSig = np.linalg.inv(Sig)    
    logPs = np.zeros(nData)
    
    const = - nDim / (2.) * np.log(2*np.pi) - 1 / (2.) * np.log(detSig) # sacalar
    for i in range(nData):
        x = X[i,:]
        expval = - 1/2. * (x - mu).reshape(1,-1).dot(invSig).dot((x - mu).reshape(-1,1)) # 2d array
        logPs[i] =  const + expval
    return logPs

Full covariance matrix

In [9]:
probC1tst1 = getLogGaussian(tstX1, estM1, estCov1)   # array
probC2tst1 = getLogGaussian(tstX1, estM2, estCov2)
probC1tst2 = getLogGaussian(tstX2, estM1, estCov1)
probC2tst2 = getLogGaussian(tstX2, estM2, estCov2)
print sum(probC1tst1 - probC2tst1 > 0)
print sum(probC2tst2 - probC1tst2 > 0)
accuracy = (sum(probC1tst1 - probC2tst1 > 0) + 
    sum(probC2tst2 - probC1tst2 > 0))/(2.*tstdatanum)
print accuracy

658
934
0.796

From graphical model

In [10]:
# Graphical model
# P(X1,...,X10) 
#        = P(X1,X2)*P(X3,...,X5|X1,X2)*
#          P(X7,...,X10)*P(X6|X3,...X5,X7,...X10)
#
#        = P(X1,...,X5)*P(X3,...,X10)/P(X3,X4,X5)
#

estM1_12345 = sum(X1[:,0:5],axis=0)/datanum
estM2_12345 = sum(X2[:,0:5],axis=0)/datanum
estM1_34567890 = sum(X1[:,2:10],axis=0)/datanum
estM2_34567890 = sum(X2[:,2:10],axis=0)/datanum
estM1_345 = sum(X1[:,2:5],axis=0)/datanum
estM2_345 = sum(X2[:,2:5],axis=0)/datanum


estCov1_12345 = np.dot(X1[:,0:5].T, X1[:,0:5])/datanum - np.outer(estM1_12345, estM1_12345)
estCov2_12345 = np.dot(X2[:,0:5].T, X2[:,0:5])/datanum - np.outer(estM2_12345, estM2_12345)
estCov1_34567890 = np.dot(X1[:,2:10].T, X1[:,2:10])/datanum - np.outer(estM1_34567890, estM1_34567890)
estCov2_34567890 = np.dot(X2[:,2:10].T, X2[:,2:10])/datanum - np.outer(estM2_34567890, estM2_34567890)
estCov1_345 = np.dot(X1[:,2:5].T, X1[:,2:5])/datanum - np.outer(estM1_345, estM1_345)
estCov2_345 = np.dot(X2[:,2:5].T, X2[:,2:5])/datanum - np.outer(estM2_345, estM2_345)


probC1tst1 = (exp(getLogGaussian(tstX1[:,0:5], estM1_12345, estCov1_12345))*
    exp(getLogGaussian(tstX1[:,2:10], estM1_34567890, estCov1_34567890))/ 
    exp(getLogGaussian(tstX1[:,2:5], estM1_345, estCov1_345)))
probC2tst1 = (exp(getLogGaussian(tstX1[:,0:5], estM2_12345, estCov2_12345))*
    exp(getLogGaussian(tstX1[:,2:10], estM2_34567890, estCov2_34567890))/
    exp(getLogGaussian(tstX1[:,2:5], estM2_345, estCov2_345)))
probC1tst2 = (exp(getLogGaussian(tstX2[:,0:5], estM1_12345, estCov1_12345))*
    exp(getLogGaussian(tstX2[:,2:10], estM1_34567890, estCov1_34567890))/
    exp(getLogGaussian(tstX2[:,2:5], estM1_345, estCov1_345)))
probC2tst2 = (exp(getLogGaussian(tstX2[:,0:5], estM2_12345, estCov2_12345))*
    exp(getLogGaussian(tstX2[:,2:10], estM2_34567890, estCov2_34567890))/
    exp(getLogGaussian(tstX2[:,2:5], estM2_345, estCov2_345)))
GraphicalModelAccuracy = (sum(probC1tst1 - probC2tst1 > 0) + 
    sum(probC2tst2 - probC1tst2 > 0))/(2.*tstdatanum)
print GraphicalModelAccuracy

0.8955
In [11]:
print accuracy
print GraphicalModelAccuracy

0.796
0.8955